Biological evaluation, molecular modeling and dynamics simulation of phenanthrenes isolated from Bletilla striata as butyrylcholinesterase inhibitors

As part of our continuous studies on natural cholinesterase inhibitors from plant kingdom, the 95% ethanol extract from tubers of Bletilla striata showed promising butyrylcholinesterase (BChE) inhibition (IC50 = 8.6 μg/mL). The extracts with different polarities (petroleum ether, ethyl acetate, n-butanol, and water) were prepared and evaluated for their inhibition of cholinesterases. The most active ethyl acetate extract was subjected to a bioassay-guided isolation and afforded twenty-two bibenzyls and phenanthrenes (1–22). All isolates were further evaluated for their BChE inhibition activity, and five phenanthrenes presented promising capacity (IC50 < 10 μM). Further kinetic studies indicated their modes of inhibition. Compounds 6, 8, and 14 were found to be mixed-type inhibitors, while compounds 10 and 12 could be classified as non-competitive inhibitors. The potential interaction mechanism of them with BChE was demonstrated by molecular docking and molecular dynamics simulation, showing that they could interact with catalytic active site and peripheral anionic site of BChE. These natural phenanthrenes provide new scaffold for the further design and optimization, with the aim to discover new selective BChE inhibitors for the treatment of AD.


Results and discussion
Chemistry. Isolation scheme of B. striata tubers was shown in Fig. S1. Two ChEs inhibition of 95% ethanol, petroleum ether, ethyl acetate, n-butanol, and water extracts of B. striata tubers were investigated and presented in Table S1. The results were expressed in IC 50 (the half-inhibition concentration μg/mL). According to the results, EtOAc extract displayed the most potent BChE inhibition and the highest selectivity in comparison to AChE (IC 50 AChE = 224.0 μg/mL, IC 50 BChE = 2.3 μg/mL). Due to the potent BChE inhibition, EtOAc extract was subjected for further bioassay-guide isolation to explore its active phytochemicals. Consequently, 22 compounds (1-22, Fig. 1) were obtained by following combination of chromatographic techniques and assigned as phenanthrenes and bibenzyls by 1 H-, 13 C-, 2D NMR, HR ESI-MS spectrum and HPLC chromatogram analysis.
As previously reviewed, a fairly large number of naturally occurring phenanthrenes have been reported from higher plants, most of them were isolated from the species of the Orchidaceae family, in 49 species: in particular Dendrobium, Bulbophyllum, Eria, Maxillaria, Bletilla, Coelogyna, Cymbidium, Ephemerantha, and Epidendrum 28 . Our present finding concerning the chemical constituents goes in line with previous reports on B. striata. Noticeably, compounds 18 and 20 are two new phenanthrenes 29 , while compounds 17, 19 and 21 were obtained from this species for the first time. Notably, the biological activities of most of isolates have not been reported.
Biological activity evaluation. For all isolated compounds, inhibitory activities on two ChEs in comparison to the references galanthamine and tacrine were detected. The initial tests were achieved at a concentration of 25 μg/mL (final concentration in the reaction system) and all of them exhibited no or weak AChE inhibitory activity. Interestingly, BChE was found to be much more sensitive than the AChE towards the evaluation (Table 1). Subsequently, IC 50 values towards BChE were determined. Sixteen compounds displayed moderate inhibition with IC 50 values lower than 100 μM, and most of them (except 3 and 11) demonstrated better inhibition than the standard galanthamine. More importantly, five phenanthrenes (6, 8, 10, 12, and 14) displayed significant BChE inhibitory effect (IC 50 < 10 μM) and their selectivity towards BChE was superior to that of tacrine.
Notably, phenanthrenes 6-13 were found to inhibit BChE with IC 50 values ranging from 2.1 to 44.6 μM. By comparing their structure and activity, it should be noted that the presence of p-hydroxybenzyl attenuated the inhibitory capacity. Compound 10, bearing two hydroxyls and one methoxy, along with the p-hydroxybenzyl substituent at C1 position on phenanthrene skeleton, presented the strongest inhibition (IC 50 = 2.1 μM) and the highest selectivity towards BChE (SI > 34.7) among the compounds with similar structures. In addition, among the five phenanthrenes dimers, 14-16 displayed attractive activity, while 20 exerted weak inhibition potency and 22 presented moderate efficacy, indicating the influence of the linkage of two monomers and the substitution position of hydroxy/methoxy on BChE inhibition. In contrast, two phenanthrofurans (18 and 19) demonstrated obvious difference potencies, and it is speculated that the substitution sites of methoxy affect whether small molecules can penetrate deep to the bottom of the active pocket and bind to the residues of key amino acids. It is worth noting that all isolated bibenzyls showed poor anti-BChE capacity compared to phenanthrenes. Collectively, the above results demonstrated that phenanthrenes should be the main BChE inhibitory components of B. striata.

Kinetic study of selective BChE inhibitors.
To study mechanism of BChE inhibition, five most potent inhibitors were selected to further study the kinetics using Lineweaver-Burk plots. As shown in Fig. 2, with the increasing concentrations of compounds 6, 8 and 14, the V max were decreased and K m were increased, these results suggested that these three compounds were the mixed-type inhibitors where they can bind to the active site of BChE. Furthermore, with the increasing concentrations of 10 and 12, the V max were decreased and K m

Molecular docking study.
To further reveal the interaction modes of five inhibitors, the molecular docking was implemented using the Glide module of Schrödinger. The active site of BChE consists of a peripheral anion site (Asp70), a choline binding pocket (Trp82), an acyl-binding pocket (Val288, Leu286, Trp231) and catalytic residues (Ser198, His438 and Glu325) 30 . Trp82 and Ser198, in particular, are considered as key amino acid residues 31 . As depicted in Fig. 3C, compound 10 bind to the active pocket of BChE via multiple hydrogen bonds with His438, Pro285, Gly115 and π-π stacking with Tyr332 and Trp82. This finding is in good agreement with the results of kinetic study and give a reasonable explanation for its highest potency against BChE. In addition, 10 and the other four inhibitors shared similarities in binding to the enzyme, with the hydroxyl groups at the C-2 and C-7 positions (except for the C-2 hydroxyl group of 8) establishing H-bonds with residues in the active pocket of the enzyme (Fig. 3). Accordingly, compound 6 was less active owing to methoxy substitution at C-7. As mentioned above, the hydroxyl groups at C-2 and C-7 of the phenanthrene skeleton are essential for inhibitory activity. Interestingly, the p-hydroxybenzyl moiety of four ligands (6, 8, 10 and 12) engaged in hydrogen bond and π-π stacking interactions with Trp82 and His438 (Fig. 3). These observations indicated that the substituents with the donor or receptor of hydrogen bond on the C-1 and C-8 play a key role in improving binding affinity of phenanthrenes. Taken together, these docking results may explain why five phenanthrenes display the highest affinity for BChE in all phenanthrenes and bibenzyls.

Molecular dynamic simulations.
Based on these promising results obtained from molecule docking, further insights concerning the dynamic behavior and binding mechanistic information were demonstrated by molecular dynamics simulation (MD) in GROMACS. The docking pose of five inhibitors with BChE (PDB:4TPK) for 50 ns were collected to generate RMSD (Root mean squared deviation) and RMSF (Root mean squared fluctuation) graphs (Fig. 4). www.nature.com/scientificreports/ RMSD of five inhibitors were plotted to verify molecular interactions and structural stability. As depicted in Fig. 4A, all five complexes possessed stable RMSD values which were all less than 0.35 nm, indicating that the structure of each complex was relatively stable. The RMSD trajectory of BChE-6 reached equilibrium within approximately 30 ns and remained until the end of the simulation. While the BChE-8 system performed best with its RMSD rising gradually before 20 ns and stabilizing in subsequent simulations without significant fluctuations, stabilizing around 0.23 nm. Also, BChE-10 system showed a gradual increase in RMSD from the beginning, a decreasing trend at about 15 ns, and was essentially stable after 30 ns. Meanwhile, the BChE-12 system had the largest fluctuations at the beginning, but it started to stabilize around 20 ns until the end of the simulation. In particular, the BChE-14 system exhibited an increasing trend at the beginning of the simulation, then largely stabilized from approximately 20 ns onwards and fluctuated in the 0.25-0.30 nm range. In a word, the RMSD plots indicated that the five inhibitors remained stable during the MD simulation by forming molecular interactions with surrounding residues.
RMSF can be used to calculate the fluctuation amplitude of a single residue and measure the degree of freedom of atomic motion in the process of dynamic simulation, and thus reflecting the flexibility of protein region motion. The RMSF of the C-alpha for the complexes of each inhibitor with BChE were calculated and the results obtained by calculating RMSF of the five systems in the equilibrium state are shown in Fig. 4B. The average fluctuation range was in the range of 0.1-0.6 nm. In short, no significant conformational changes in protein and molecular structure were observed during simulations and the vital residues including Trp82 and Phe329 in them were quite stable. The average RMSF values obtained for the active molecules also indicated good protein flexibility. Accordingly, all five inhibitors exhibited a good stability and inhibitory activity.
The radius of gyration (Rg) correlated with the relationship between the atomic mass and the center of gravity of a protein, which can be used to evaluate the structural compactness of the protein docking complexes. The stably folded proteins showed a relatively steady value of Rg, while the perturbation of protein conformation showed fluctuations in Rg values. As shown in the Fig. 5A, the Rg value of both (10 and 12) remain stable in the range of 2.30-2.35 nm. From the relative frequency analysis (Fig. 5B), the Rg values of protein complexes were slightly higher than those of individual BChE. But the offset value of Rg is only 0.01-0.015 nm. This suggested that the protein BChE became marginally looser bound to the compounds 10 and 12, but the effect is limited.
The solvent-accessible surface area (SASA), an important parameter to describe the hydrophobicity of proteins, was observed and shown in Fig. 6A. Results showed that SASA value of both (10 and 12) remain relatively stable in the range from 220 to 245 nm 2 in the simulation time 0-50 ns. The relative frequency analysis (Fig. 6B) In silico ADME prediction. To evaluate drug-like properties of compounds, the ADME (absorption, distribution, metabolism and excretion) properties of five inhibitors were predicted in silico using the QikProp v. 5.5 (Schrödinger). As showed in  were shown in sticks colored by atom type, the carbon in green, hydrogen in white, nitrogen in blue, and oxygen in red. Interactions included non-covalent bonds and π interaction, the hydrogen bonds in yellow, π-π stacking in blue. www.nature.com/scientificreports/ solubility of the molecules. The QPlogPo/w values of five ligands were lower than 5.311, signifying that they had good water solubility. Rather, the QPlogBB values which predict the ability of molecules to cross the blood-brain barrier were in an ideal range, indicating that five molecules could enter the central nervous system via oral administration. In addition, human oral absorption (HOA%) values exceed 90%, implying that these molecules were fully absorbed in the gastrointestinal tract.
After the determination of spontaneous hydrolysis rate of ATCI/BTCI, 25 μL of 0.226 U/mL ChE (AChE from electric eel and BChE from equine serum) in PBS (pH 8.0) were added to each well of the 96-well plates and then the plates were shaken 15 s under the medium speed. Finally, changes in the absorbance at 405 nm were determined by the method mentioned above to obtain the apparent enzymatic hydrolysis rate. The actual enzymatic hydrolysis rate was corrected by deducting the spontaneous hydrolysis rate of ATCI/BTCI. The enzymatic hydrolysis rate determined when 10 μL of methanol added to the well instead of the samples were set as blank. The concentration of the methanol in final reaction mixture is lower than 4% and has no influence on AChE/ BChE. All the tests were achieved in triplicate and compounds showing more than 50% inhibition rate at 25 μg/ mL (final concentration in the reaction system) were further determined the IC 50 values. The corresponding IC 50 measurement concentration were set based on the inhibition rate at 25ug concentration, ranging from 0.04 to 50 ug/ml. The least square method of inhibitor versus Response-variable slope (four parameters) model in GraphPad Prism 8.0 (GraphPad Software, USA) was used for calculation and fitting the IC50 values. The inhibition percentage of ChE (%) calculated using the following formula: V blank and V sample represent the hydrolysis rate of the blank control and test sample respectively.

Kinetic characterization of BChE inhibition.
The Lineweaver-Burk plots analyses was performed to determine the inhibition mode of five potential inhibitors (IC 50 < 10 μM). Five increasing concentrations of BTCI as a substrate in the presence of three different concentrations of compound were used in the kinetic measurement. The concentrations of the enzyme and DNTB used in the determination are 0.226 U/mL and 3 mM respectively. The inhibition constants (K i ) of the inhibitors were determined by the secondary plots which constructed using the slopes of the Lineweaver-Burk plots.
Docking study. The molecular docking of active compounds was performed using the Glide module of Schrödinger. The X-ray crystal structure of human BChE (PDB ID: 4TPK) was obtained from Protein Data Bank. The five compounds were prepared by LigPrep module. The Protein Preparation Wizard of Maestro was used to prepared the complex structure of BChE, which retain one protein conformation, supplement incomplete residues, add hydrogen and keep the active site ligand only. The pH value of Protein and Compounds were set as the default values 7.0 ± 2.0. Then we used the Receptor Grid Generation module to define the active site of the protein through the cocrystallized ligand, and it was excluded from the grid generation. The redocking of cognate ligand was performed to evaluate the reliability of the docking model and the RMSD value was 1.95 Å (< 2 Å), which means the docking model is credible. The other parameters were set at the default values in Schrödinger.

Molecular dynamic simulations.
To further understand the mode of action of the inhibitors in the biological environment as well as to test the stability of the complexes with the hit compounds, MD simulation was carried out on GROMACS 2019.5. It was employed to obtain the trajectories of five complexes, which were used to obtain root mean square deviations and fluctuations (RMSD/RMSF) thus to evaluate the stability of each complex. The molecular dynamics experiment of 10 and 12 were repeated three times.
The topology of the protein was generated by GROMACS using the GROMOS96 force field, and the topologies of the ligands were obtained from the Automated Topology Builder website at https:// atb. uq. edu. au/ index. py. To maintain the periodic boundary conditions, the minimum distance between the protein atoms and the cell wall was set to 1 Å throughout the MD simulation process with the cubic cell. The simple point charge solvent model (spc216) was employed to solvate the system, and the species of the solvent was water. The processed system might not be in the electrically neutral state, and the charge of the system could be balanced by adding counter-ions such as Cl − or Na + to achieve the physiologically neutral state. The energy of the system was minimized to eliminate steric clashes, create an index and modify the .mdp file for the thermal coupling group as well as the subsequent position limitation. The position restraint dynamic under NVT (constant volume) and NPT (constant pressure) conditions at 300 K were brought into operation by a V-rescale which was modified by a Berendsen thermostat. The Particle Mesh Ewald method was employed to compute the long-range Coulombic www.nature.com/scientificreports/ and Lennard-Jones interaction energies for the protein-ligand interactions and the ligand dynamic. Different modules in the GROMACS package were brought into effect to analyze the results of the MD simulation to acquire the RMSD and RMSF. The simulation method used in this study is the same as that previously described in our previously published study 49 . ADME prediction. ADME properties are all used to assess the drug-forming properties of drugs. QikProp is a module in Schrödinger for predicting ADME properties of molecules. QikProp contains two basic models, fast and normal. The results were obtained for a variety of properties including partition coefficients (QPlogP octanol/water), predicted water solubility (QPlogS), percentage oral absorption in humans, volume, hydrogen bond donors and acceptors, etc. Forty-four properties of the molecule can be obtained after prediction. Before performing the QikProp calculations, molecular pre-processing was completed and the parameters were set as default.
Statistical analysis. All analyses were carried out in triplicates. Data obtained were presented as mean and standard deviation values (S.E.M.). The IC 50 values were calculated with GraphPad Prism 8.0.

Data availability
All data is provided in the manuscript or supplementary material.